Effect of confinement on dense packings of rigid frictionless spheres and polyhedra 
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We study numerically the influence of confinement on the solid fraction and on the structure of 
three-dimensional random close packed (RCP) granular materials subject to gravity. The effects 
of grain shape (spherical or polyhedral), material polydispersity and confining wall friction on this 
dependence are investigated. In agreement with a simple geometrical model, the solid fraction 
is found to decrease linearly for increasing confinement no matter the grain shape. Furthermore, 
this decrease remains valid for bidisperse sphere packings although the gradient seems to reduce 
significantly when the proportion of small particles reaches 40% by volume. The confinement effect 
on the coordination number is also captured by an extension of the aforementioned model. 



I. INTRODUCTION 

Granular materials are well known for their wide range 
of fascinating properties. Their theoretical description is 
difficult for many reasons. One of them is the importance 
of the local arrangement of grains within the material on 
its macroscopic behavior. Real granular systems have 
boundaries, but for the sake of simplicity, scientists often 
neglect them (see for example 0, Q ) and the system is 
then considered as infinite. This assumption is not al- 
ways justified since the presence of boundaries modifies 
the system local arrangement in their vicinity. Moreover, 
due to the intrinsic steric hindrance of granular mate- 
rials those structure modifications often propagate over 
distances in the order of several grain sizes. As a conse- 
quence, the behavior of granular systems may be strongly 
influenced by the presence of sidewalls even if the con- 
finement length is large compared to the grain size. 
The crucial role of confinement on system properties 
has been pointed out in many works dealing with grav- 
ity driven granular flows [3|-|8|, granular segregation l9[ , 
structure and mechanics of granular packings [10l - [l2l ). 
granular systems in narrow silos (l~3| or granular pene- 
tration by impact [l4j|. Those studies point out that two 
major physical properties can be influenced by the pres- 
ence of sidewalls. First, they can induce friction that 
might be important in respect to the internal friction of 
the system Q explaining the well-know Janssen effect [l5[ 
or unexpectedly hig h angle values observed with confined 
granular heaps jl6l [l7j or confined chute flows 0]. Sec- 
ond, as mentioned above, they might also alter the ge- 
ometrical structure of the system near the wall, where 
particles tend to form layers, giving rise to a fluctuat- 
ing local solid fraction with distance from the wall [TlT ] 
and affecting the properties of confined systems. Note 
that the effect of confinement is not limited to the vicin- 
ity of the walls but may propagate within the sample. 
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This is more particularly the case for confined granular 
chute flows for which it has been shown that the good 
dimensionless number to quantify the sidewall effect is 
not the number of grains per unit of width between side- 
walls but the ratio of the flow height to the gap between 
sidewalls 

Here, we focus mainly on the geometric effect of the 
presence of sidewalls on quantities like the solid fraction 
and the coordination numbers. Recently, Desmond and 
Weeks used numerical simulations to study the effect of 
confinement on binary atomic systems at the random- 
close-packing limit [ 1 81 ] - Their numerical results agree 
with a simple geometrical model fl9U2l| (hereafter called 
the geometrical model) which captures the evolution of 
the solid fraction of random close packings of spheres 
with confinement. It is based on the following config- 
uration: a packing of particles is confined between two 
parallel and flat walls separated by a gap W . It then as- 
sumes that such a confined system is made of two bound- 
ary layers (of thickness h) and a bulk region and that the 
solid fraction of the boundary layers, 4>bl, is lower than 
that of the bulk region <pbuik- By writing the total solid 
fraction <\> as the average of both the bulk region and 
boundary layers solid fractions (weighted by their rela- 
tive thickness), the geometrical model predicts that the 
average solid fraction decreases linearly with 1/W: 
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where C — 2h{(j>b u ik — 4*bl)- Note that this model can 
be easily adapted to other boundaries such as cylindrical 
ones [18[ . The three parameters of the geometrical model 
{4>buik, 4>bl and h) probably depend on grain shape, 
packing polydispersity and confining wall properties. 
Here we study the effect of confinement on quasi-static 
dense frictionless granular systems (i.e. grains interact- 
ing through hard core repulsion) subjected to gravity. 
We test the validity of the geometrical model for such 
systems and study the aforementioned dependencies. 
Using numerical simulations, we investigate the actual 
effect of grain shape by comparing packings of spheres 
with packings of polyhedra. Furthermore, we assess the 
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effect of packing polydispersity by comparing monosizcd 
and binary packings. Besides, we check the effect of 
grain-wall friction. Eventually, we look into packing 
microstructure by studying the effect of confinement on 
the coordination number. 

The paper is organized as follows. Section HT1 describes 
our simulations with details as well as the numerical sim- 
ulation method used. After checking the state of pack- 
ings in section ITO1 section HVl investigates how the solid 
fraction is modified by confinement and how these mod- 
ifications are influenced by packing polydispersity, grain 
shape and confining wall friction. Then, we examine in 
section [V] the modification of the packing microstructure 
with confinement. Finally, in section l*VTl we summarize 
our results and present our conclusions. 

II. SIMULATION METHODOLOGY 

A. Geometry of grains 

The simulated system is a three-dimensional dense as- 
sembly of n frictionless rigid grains of mass density p, 
interacting with each other through totally inelastic col- 
lisions. 

Since grain shape may influence the behavior of granu- 
lar materials [22U25| , two types of grains have been stud- 
ied: spherical grains of average diameter d and polyhc- 
dra of average characteristic dimension d. The polyhcdra 
shape (Fig. [T]) is that of a pinacoid, with eight vertices, 
fourteen edges and eight faces. This polyhedron has three 
symmetry planes and is determined by four parameters: 
length L, width G, height E and angle a. According 
to an extensive ex per imental study with various rock 
types reported by [2g, the pinacoid gives the best fit 
among simple geometries for an aggregate grain. In or- 
der to have similar aspect ratio for both grain shapes, 
the pinacoid dimensional parameters were taken identi- 
cal (L = G = E), with the characteristic dimension d 
expressed as d = L. In addition, angle a was set to 60°. 
For each grain shape, two grain diameters (or charac- 
teristic dimensions) have been considered: large g?l and 
small dg = eft/2. 

B. Samples preparation 

The packing geometry is that of a parallelepiped 
(Fig. (2) of dimensions L x by L y by L z . Periodic Bound- 
ary Conditions (PBC) are applied in the x direction to 
simulate an infinitely long parallelepiped using a finite 
number of grains. The packing is confined in the y di- 
rection between two fixed parallel walls separated from 
each other by a L y = W large gap. In some cases, PBC 
are also applied along the y axis to simulate unconfincd 
reference state, with W set to 20cft. The packing is sup- 
ported on the xy-plane by a fixed frictionless bottom wall 
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FIG. 2. (Color online) Typical 3D snapshots of packings made 
of polysized spheres (a) and pinacoids (b). The confinement is 
characterized by the gap W between sidewalls. The direction 
of gravity is —z. 
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and delimited by a free surface at its top. 

Grain samples are composed of various proportions of 
small and large grains having the same shape. In order to 
reduce the thickness of the crystallized layer commonly 
observed inside confined packings at the interface with 
smooth walls [13, HH, each population of grain size is 
randomly generated with a Gaussian distribution char- 
acterized by its mean d and its variance d 2 /900. Any- 
how, for the sake of simplicity, packings made of a unique 
population of grains (either small or large) will be called 
" monodisperse" , whereas packings made of small and 
large grains will be called "bidisperse" in the following. 
In the latter case, the proportions of small (xs) and large 
(xl) grains expressed as percentages by volume are of 
course linked through xs = 100 — xl- 

Each sample is constructed layer- by- layer according to 
the following geometrical deposition protocol inspired by 
[29l ]: spherical particles of a sufficient diameter to cir- 
cumscribe the larger grains of the sample are sequen- 
tially dropped along z in a parallelepiped box having the 
aforementioned geometry. Each particle stops on the free 
surface made of the underneath layer of particles (or on 
the bottom wall for the first layer of particles) and is fur- 
ther moved so that it lays on three particles chosen to 
locally minimize its altitude z. Finally, the sample ac- 
tual grains are randomly substituted for those spherical 
particles. For polyhedra samples, a random orientation 
is further assigned to each pinacoid. Note that according 
to this protocol, some of the deposited grains may not be 
in contact with their neighbors depending on their size 
and shape. 



C. Initialization and solicitations 

The system initialization is identical for spheres and 
pinacoids samples. The first step consists in geometri- 
cally depositing n grains into a parallelepiped box, then 
PBC are substituted for the lateral walls of the paral- 
lelepiped box along the x direction (along the x and y 
directions for biperiodic reference state). Finally, gravity 
g (0,0,— g) is applied in order to compact the sample. 



normal and tangential relative velocities at contact point 
respectively before and after collision, the formal normal 
and tangential relative velocities are defined as follows: 
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where ejv and ej~ measure the inelasticity of collisions 
and reduce to the normal and tangential restitution co- 
efficients in the case of binary collisions. 

These generalized contact laws support momentum 
propagation through contact networks inherent to dense 
assemblies of grains. For a given time step, impulses and 
velocities are determined according to an iterative pro- 
cess using a non- linear Gauss-Seidel like method [33 |. 
In the case of large size packings of rigid grains, the CD 
method supports larger time steps, leading potentially to 
faster calculations than the molecular dynamic method 
for which small time steps are needed. 

The CD method was applied using the LMGC90© 
platform [33l . HiJ which namely implements a 3D ex- 
tension of a 2D contact detection algorithm described 
with details in [35| . Basically, contacts with a given 
grain are sought exclusively among its neighbors. When 
a neighbor is located closer to the grain than a thresh- 
old distance called gap, a 3D extension of the shadow- 
overlap method devised by Moreau [33l . |3(| is applied. 
In case of overlap between the grain shades, their con- 
tact plane is determined. Four contact situations may 
be encountered (Fig. [3]): vertex-to-face or edge-to-edgc, 
represented by a single point and called simple; edge-to- 
face, represented by two points and called double; finally, 
face-to-face, represented by three points and called triple 
(vertex-to-edge and vertex-to-vertex being very unlikely 
to happen). These situations allow identifying a con- 
tact plane and compute the contact impulse and velocity 
components at each contact point. 

This method proved apt to deal with dense flows of 

M. 
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tion [35l . 



as well as with quasi-static plastic deforma- 
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D. Contact dynamics method 

Discrete numerical simulations were performed using 
the contact dynamics (CD) method [13, HU, which is spe- 
cially convenient for rigid grains. This method is based 
on implicit time integration of the equations of motion 
with respect to generalized non-smooth contact laws de- 
scribing non-interpenetration and dry friction between 
grains. This formulation unifies the description of last- 
ing contacts and collisions through the concept of im- 
pulse, which can be defined as the time integral of a force. 
The generalized non-smooth contact laws are exp ress ed 

in terms of impulse Pc and formal relative velocity Vc at 
contact point C. IfV, 
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FIG. 3. (Color online) Two polyhedra can experience simple 
contacts (a), double contacts (b) or triple contacts (c). 



E. Materials and system parameters 

The present study focuses on monodisperse sphere 
packings (MSP), bidisperse sphere packings (BSP) and 
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monodisperse pinacoid packings (MPP). 

The spacing of lateral walls W takes discrete values 
between 5cIl and 20d_L, and the sample period along the 
x axis is 20g?l- With the final height h of the packing in 
the range of 16di to 20g?l, the number n of grains varies 
between 1,900 and 30,400 for spheres, depending on the 
proportion by volume of small grains, and between 3,600 
and 15,000 for pinacoids. 

The time step value At is taken small enough to mod- 
erate the grain interpenetration incumbent to grains dis- 
placement between two successive implementations of the 
contact detection algorithm, but sufficiently large to keep 
the calculation duration reasonable. In this perspective, 
limiting to c^/lOO the translation of grains during At 
seems appropriate. For our grain packings subject to 
compaction under their own weight, this leads to the fol- 
lowing relation: 

= Vmax- At, (3) 

with v max = \J1gAh the maximum speed reached by 
a grain free falling from initial height of the deposited 
packing down to the altitude of the packing free surface 
at the end of the compaction. Hence, equation [3] leads to 
the following expression for the time step: 

m= j_. x [aL.K. (4) 

100 V 2Ah y g { ' 

Although the between-grain friction is set to zero, that 
of wall-grain contacts (/im) is assigned non-zero values in 
a few simulations to study the influence of wall friction. 

The simulated system parameters are summarized in 
table U for spheres and in table [TT] for pinacoids. They 
are expressed as dimensionless quantities by defining the 
following normalization terms: lengths and times are re- 
spectively measured in units of and to = cLl /g, the 
characteristic free fall time of a rigid grain of diameter 
subject to gravity exclusively. For a given set of system 
parameters, three grain packings are simulated in order 
to average the various measured quantities. 



TABLE I. Sphere packings 
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TABLE II. Pinacoid packings 
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packings for various grain shape and polydispersity, it 
is necessary to adopt a reference packing state and to en- 
sure that the compaction method used allows to approx- 
imate such a state while providing sufficient repeatabil- 
ity for a given set of materials and system parameters. 
As mentioned in paragraph lll CI the compaction method 
used consists in depositing rigid frictionlcss grains (with 
or without wall friction) under their own weight. For 
sphere packings with presumably negligible confinement, 
several authors have experimentally [45l [4(| or numeri- 
cally 0, H(| HtJ observed that this compaction method 
led to random close-packed states characterized by the 
generally agreed solid fraction value of 0.64. Accord- 
ing to j48[, random close- packed states of rigid friction- 
less grains (spherical or non-spherical) are equivalent to 
packing states in which the grains are homogeneously 
spread and in a stable equilibrium without crystalliza- 
tion or segregation (observe that the notion of "stable 
equilibrium" refers to the minimization of a potential en- 
ergy that ensures maximum solid fraction [49j| ) . Further- 
more, extensive investigation of the random close-packed 
state carried out by [2J with spherical particles has ev- 
idenced the uniqueness of this state in the limit of in- 
finitely large samples subject to fast isotropic compres- 
sion (to avoid cristallization). Hence, the influence of 
wall-induced confinement on the solid fraction and struc- 
ture of dense packings may be assessed against the ran- 
dom close-packed state taken as the reference. Keeping in 
mind that our compaction method only allows to approx- 
imate the random close-packed state (our compression is 
not isotropic) and that the uniqueness of this reference 
state has only been evidenced for sphere packings, it is 
expected that meeting as much as possible the criteria 
stated by [48[ will lead to sufficiently repeatable solid 
fraction and microstructure characteristics for a given set 
of materials and system parameters to observe confine- 
ment effects for various grain shapes and polydispersity. 
Therefore, preliminary assessment consists in checking 
the state of our simulated packings (both sphere and 
pinacoid packings) in terms of stable equilibrium, homo- 
geneity and reasonable interpenetration given the partic- 
ularities of the contact dynamics method. Further as- 
sessment will be undertaken in sections [IV] and [V] 



III. STATE OF PACKINGS 



A. Equilibrium 



In order to examine the influence of wall-induced con- 
finement on the solid fraction and structure of dense 



According to Q, grain packings for which the follow- 
ing criteria are met on each grain have reached a stable 
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equilibrium: 



F < ICtVP, 



J2 M < 10- 4 d 3 P, 



E r < 10~ 8 d 3 P. 



(5) 



(6) 



(7) 



where E F, E M and E c are respectively the net force, 
net momentum and total kinetic energy of the grain. 
Indeed, Q have observed that setting to zero all grain 
velocities in such a state and letting the packing relax 
further did not lead to any kinetic energy or unbalanced 
force level regain beyond these threshold values. 

As a consequence, these criteria were used to check 
the attainment of a stable equilibrium state in our simu- 
lations, which was the case for all of them. 



B. Interpenetration 

The grain interpenetration, calculated as the sum of 
interpenetrated volumes between neighboring grains di- 
vided by the sum of grain volumes, was checked in the 
bulk region of the packing at the end of each simulation. 

For sphere packings, the interpenetration was calcu- 
lated analytically as the sum of interpenetrated volumes 
between couples of spheres (for a given couple of spheres, 
the interpenetrated zone consists of two spherical caps) 
and it was found to be very low (in the range of 1CD 5 to 
1(D 3 % by volume). 

For pinacoids, a routine was designed to compute the 
solid fraction as well as lower and upper bounds of the 
grain interpenetration. Basically, this routine consists of 
superimposing a lattice on the grain packing and calcu- 
lating the solid volume in each cell of the lattice. For 
a given cell, this solid volume is the sum of elementary 
volumes analytically calculated from the intersection be- 
tween any pinacoid and the cell. In order to bound the 
grain interpenetration, one shall focus on cells intersected 
by two neighboring pinacoids, leading to one of the four 
situations depicted on Fig. @] (in 2D for simplicity rea- 
sons): 

• In situations a and b, the solid volume V so i con- 
tained by the cell is in excess of actual cell vol- 
ume Vceii, hence the lower bound of actual inter- 
penetrated volume is (V so i — V ce u)/2 (situation a) 
whereas the upper bound is V so i/2 (situation 6); 

• In situations c and d, the solid volume contained 
by the cell is smaller than actual cell volume, hence 
the lower bound of actual interpenetrated volume 
is (situation d) whereas the upper bound is V so i/2 
(situation c). 
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FIG. 4. 2D simplified representation of the four intersection 
situations between two pinacoids and a lattice cell: (a) solid 
volume > cell volume and partial interpenetration in the cell; 
(b) solid volume > cell volume and total interpenetration in 
the cell; (c) solid volume < cell volume and interpenetration 
(partial as represented or total); (d) solid volume < cell vol- 
ume and no interpenetration. 



Observe that V so i/2 is the upper bound of the inter- 
penetrated volume no matter the situation. When their 
size decrease, the lattice cells that are intersected by 
two pinacoids tend to concentrate exclusively inside 
actual interpenetrated areas (I A) where V so i/2 — V ce u 
or astride their border (AB) where < V so i/2 < V ce u. 
Hence, the total interpenetrated volume of the packing 
Vi is bounded by the following interval: 



Vj € [ E Vcell 
cellelA 



cellelA 



Yl ^cell] 
elleAB 



in which E Vcell tends to zero with decreasing cell 

cellGAB 

size. 

For each pinacoid packing geometry, table Mil gathers 
lower I-min and upper I max bounds of grain interpenetra- 
tion, e.g bounds of Vj expressed as a percentage of the 
packing solid volume. These values were computed in 
the bulk region using a lattice with d^ /20-large cubical 
cells, and each of them was averaged over three simu- 
lations. The interpenetration calculated in our pinacoid 
packings, in the range of 3 to 5% by volume, is clearly 
much higher than the one calculated for sphere packings. 

A first reason to explain these differences lies with the 
determination of contact between two grains. In the case 
of sphere packings, this determination is very simple and 
requires no interpenetration: grains are in contact when 
the distance between their centers is lower or equal to 
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TABLE III. Pinacoid packings interpenetration for MPP. 
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3.7 



the sum of their radii. Such a contact is only one point, 
which is located on the segment connecting the centers 
of spheres at a distance of each sphere center equal to 
its radius. Besides, the orientation of the contact nor- 
mal is borne by the segment connecting the grain cen- 
ters. In the case of pinacoid packings, the determination 
of contact between two grains is much more complex, 
time-consuming and implies more or less interpenetra- 
tion: first, grains are in contact when their respective 
shadows always overlap no matter the projection direc- 
tion. Hence, much more calculation than for sphere pack- 
ings shall be performed to prove the existence of a con- 
tact, and the simultaneous achievement of these overlap 
situations generally implies some interpenetration. Next, 
in case of a contact, it may not be a unique point but 
rather two (edge-to-face contact) or three (face-to-face 
contact) points as explained in section l"HD| which obvi- 
ously leads to more interpenetration. 

A second reason lies with the non-smooth approach 
associated with the contact dynamics. In molecular dy- 
namics (Hoj . contact forces increase proportionally to a 
power function of the interpenetration, leading to high 
repulsion contact forces thus low interpenetration in the 
limit of rigid grains. In the contact dynamics method 
where no such relation is applicable, the interpenetration 
is namely monitored by the quality of the convergence of 
impulses and velocities at contact points within the range 
of permissible values on the generalized non-smooth con- 
tact laws. Hence, in addition to an appropriate time step 
value, a low level of interpenetration requires optimizing 
both the convergence criteria and number of Gauss-Seidcl 
iterations while keeping the calculation time acceptable 
(for more information, refer to (30l. l3lL l33l. l34j ) . 

Anyhow, the contact dynamics method is known to 
give interpenetration values in the order of a few percent 
by volume (see HH), and our quest of the densest possi- 
ble disordered packing made of frictionless particles un- 
surprisingly leads to interpenetration values in the higher 
range. Hence the interpenetration evidenced by our re- 
sults is acceptable. 



C. Homogeneity of distribution 

In order to ensure that the applied compaction method 
leads to homogeneously distributed packings, we examine 
the variations in the proportions of large {Pl) and small 
(Ps) grains with distance z from the bottom wall (Fig. [3] 
(a) and (b)). Therefore we count the number of parti- 
cles of each size in d^-thick regions of the packing and 



divide that number by the total number of grains. Al- 
though small deviations (that may be due to segregation) 
close to the bottom of our packings are observed, the pro- 
portion profiles are almost constant showing that grains 
in sphere or pinacoid packings are reasonably vertically 
homogeneously distributed. The absence of segregation 
along the y axis is also checked for BSP in the homo- 
geneous zone (e.g. far from bottom and free surface). 
Fig. [5] (c) displays the variations of proportions Pl and 
P s in G?£-thick layers parallel to the sidewalk. The pro- 
portion profiles are almost constant, showing reasonable 
horizontal homogeneity. 



IV. SOLID FRACTION 

A. Average solid fraction 

In this subsection, our aim is to study the effect of con- 
finement on the solid fraction of MSP, BSP and MPP, 
that is to say for various proportions of small particles 
and various grain shapes. For this purpose, we report 
the evolution of the aforementioned quantity for several 
values of gap between sidewalls. We will also test the geo- 
metrical model mentioned in the introduction (cf Eq. [1} . 
The solid fraction is computed from analytical calcula- 
tion of the volume of each sphere or each pinacoid present 
within a given volume. This volume incorporates any 
particle located 3di away from the bottom wall and the 
free surface. For the solid fraction o f sp here packings, 
the use of the Voronoi tessellation [EU, [Hf gives the same 
results. 

Figure [5] reports the average solid fraction for BSP, 
MSP and MPP versus d^/W. A first observation is that 
for a fixed d^/W value, an addition of small grains in 
a monosized sphere packing increases the solid fraction. 
This well-known phenomenon can easily be understood 
by considering two limit cases. The first one consists 
of a monosized sphere packing to which we add a few 
small particles (x s <C 100%). In this case, small grains 
partially fill the porosity of the monosized packing and 
increase the solid fraction. The second limit case corre- 
sponds to a packing of small grains to which we add a 
few large particles (x s « 100%). The largest particles 
can then be considered as islands in a sea of small grains 
whose solid fraction is equal to that of a monosized pack- 
ing: 4> m ono- Since the solid fraction of the islands is equal 
to 1, the average packing fraction is greater than </> mono - 
More interestingly, an excellent agreement between our 
data and the geometrical model is found. The corre- 
sponding values of C and (fibuik are reported in tabic IIVI 
It should be pointed out that the value of <f>buik obtained 
for MSP is consistent with that of the random close pack- 
ing (0.64) reported in the literature p34| . 

Note that in Ref. [l8|, Desmond and Weeks compare 
the geometrical model with simulations of bidispcrsc 
sphere packing (50-50 binary mixture with particle size 
ratio of 1.4) in the absence of gravity. Our results show 
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FIG. 6. (Color online) Plot of the average solid fraction versus 
d L /W for MSP, BSP and MPP. The lines are fits from the 
geometrical model [Eq. [1]. Error bars denote the standard 
deviation. 
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TABLE IV. Values of C and 4>buik obtained by fitting the data 
reported in Fig. [S]with equation JT]). 
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FIG. 5. (Color online) (a) Homogeneity of MSP (o) and MPP 
(x) along the z axis, expressed as the ratio of the number of 
grains to the total number of grains in di-thick layers, (b) 
Homogeneity of BSP along the z axis, expressed as the ratio of 
the number of large (empty symbols) and small (full symbols) 
grains to the total number of grains in di-thick layers for x s = 
10% (□), x s = 25% (o) and x s = 40% (A), (c) Homogeneity 
of BSP along the y axis, same calculation method, same key. 
For (a), (b) and (c), the gap between sidewalls is W = 5dz, 
and the data have been averaged over three simulations. 



that the validity of this model is much broader since it 
still holds in the presence of gravity for monodisperse 
sphere packings, for bidisperse sphere packings (inde- 
pendently of the fraction of small grains) as well as for 
monodisperse pinacoid packings. This result is impor- 
tant in the framework of real granular materials whose 
grains are far from being perfect spheres. Let us recall 
that the fit parameter C is equal to 2h(<pbuik — 4*bl) 
(see section UJ. Our results show that when the fraction 
of small grains, x s , increases, C decreases. This can 
be the consequence of a decrease of the distance of 
propagation of the sidewall effects h or/and of the 
difference 4>Buik — 4>bl- To address this point we will 
study the local variation of the solid fraction close to the 
sidewalls. This is the objective of next subsection. 



B. Solid fraction profiles 

Figure [7] depicts the solid fraction profile as a function 
of the distance y/d,L to the left sidewall for MSP (a), 
BSP (b) and MPP (c). The local solid fraction fluctu- 
ates with the distance from the wall, especially in the 
neighborhood of sidewalls and, if W is large enough, it 
reaches a uniform value away from the sidewalls. The in- 
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FIG. 7. (Color online) Solid fraction profiles as a function 
of distance yjdh from confining wall for (a) MSP with W = 
4d L , W = 5d L , W = 10d L and W = 20d L , (b) BSP with 
W = 5d L and x a = 0%,10%,25% and 40%, and (c) MPP 
with W = 5d L , W = lOdz, and W = 2Qd L - Fluctuations of 
the local solid fraction are due to the layering of particles in 
the vicinity of the sidewalls. The inset in Fig. [TJa) is a zoom 
over 3di . The inset in Fig. [TJb) shows the solid fraction for 
BSP (here x a = 40%) and the corresponding fit (Eq. El). 



set of Fig. UJa) reports the packing fraction fluctuations 
as a function of the non-dimensional distance from the 
wall y/dL. The aforementioned fluctuations clearly re- 
flect the layering due to the presence of sidewalls - i.e. 
an order propagation in the y direction ■ For MSP, the 
confinement effect propagates over approximately 3dt to 
4di. As a result, packings for which W < 6<1l to 8<1l are 
influenced by the presence of walls over their full width. 
In other words, for such size, the order generated by the 
sidewalls propagates in the whole packing. On the con- 
trary for BSP as well as for MPP, the propagation seems 
to be shorter (approximately 1.5di to 2d,L for BSP and 
about 2d,L for MPP). The presence of bidispersity or non- 
sphericity induces disorder in the vicinity of the sidewalls 
which mitigates the layering. To quantify more precisely 
the sidewall effects we have fitted the solid fraction pro- 
files reported in Fig. [Jj by the following empirical law: 
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(8) 



In this expression, the characteristic lengths of the side- 
wall effect propagation for large (L) and small (S) grains 
are respectively A^ and A5. Parameter a characterizes 
the solid fraction increase close to the sidewalls, 7, and (3i 
respectively correspond to the period and amplitude of 
the structuration oscillations caused to the solid fraction 
profile by the layering of small and large particles. For 
monosized packings, we use the aforementioned fit with 
fis = 0. An example of the obtained fits is plotted in 
the inset of Fig. [7b. Let us stress out that the fit used 
is purely empirical. Our aim is to obtain a reasonable 
approximation for confinement effect propagation rather 
than a precise description of the solid fraction profiles by 
an equation. 

The values of Xl (those of As are not statistically 
relevant for xs < 0.25), normalized by the average 
grain size d avg , obtained this way are reported in Fig. [5] 
for W = lOd^. For sphere packings, the normalized 
characteristic length is found to decrease when the 
fraction of small spheres x s increases. Indeed, for x s = 
we have \L/d avg ~ 1.4 whereas Xh/davg ~ 0.85 for 
x s = 40%. This decrease proves that the polidispersity 
mitigates the confinement effect. Morover, the fact 
that \L/d aV g decreases with x s demonstrates that Al 
decreases quicker than the mean grain size. For MPP 
we obtain X/d avg — 1.2 which is smaller than the value 
obtained for MSP. This indicates that the sidewall effect 
is also mitigated by an increase in grain angularity. 
Hence, characteristic length Xl is expected to correlate 
with the thickness h of the boundary layers introduced 
in equation [TJ In the inset of Fig. [8] we report Xl 
versus C= 2h(4n,uik — 4>bl) and observe a good linear 
correlation between these two parameters. Furthermore, 
the data for both sphere and pinacoid packings collapse 
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FIG. 8. (Color online) Characteristic length of confinement 
effect Xl versus fraction of small grains x s for MSP and BSP 
(circle) as well as for MPP (diamond) with W = lOd. That 
length is defined in equation [8] The effect of confinement 
is found to decrease with increasing grain polydispersity or 
grain angularity. The inset reports the same length versus 
that of the fit parameter in eq. ([T]): C = 2h(4>buik — &bl)- 
The dashed line corresponds to a linear fit. 



on the same straight line whose intercept is equal to zero. 



C. Effect of grain/ wall friction 

So far, the presented simulations were performed with 
frictionless grains and sidewalls. However, additional 
simulations were performed to investigate the influence 
of grain/wall friction. For this purpose, the friction coef- 
ficient between grains was kept equal to zero, whereas the 
grain/ wall friction coefficient /i gw was successively set to 
0.3, 0.5 and 1. As before, three grain packings were sim- 
ulated for each grain/ wall friction coefficient in order to 
average the measured quantities. Our aim is not to ad- 
dress this point, but just to mention that in our contact 
dynamic simulations, we found that the grain/wall fric- 
tion had no effect since neither the average solid fraction 
nor the solid fraction profiles were affected by [A gw . This 
result demonstrates that the influence of confinement on 
packing fraction is purely geometrical. 



V. PACKING MICROSTRUCTURE 

Section HlTl has established that our packings are homo- 
geneous and that they have reached a stable equilibrium 
with acceptable interpenetration. Then, in section [IV] we 
have verified that the simulated compaction method al- 
lows to accurately achieve the 0.64 solid fraction charac- 
teristic of the random close-packed state of monodisperse 
sphere packings when PBC are substituted for sidcwalls. 



Moreover this method is sufficiently repeatable to show 
significant influence of the confinement on the solid frac- 
tion of various grain packings. Now, section [V] focuses 
on the internal state of our packings in order to investi- 
gate the influence of confinement on their microstructurc. 
We first investigate the presence of textural order (sec- 
tion |V3J. Then, we study a usual characteristic to de- 
scribe the microstructurc of grain packings: the mean 
number of contacts per grains (coordination number). 
For various grain shapes and polydispersities, section lVBl 
discusses the influence of confinement on that character- 
istic. 



A. Order 

In this section, our aim is to investigate the presence 
of long-range textural order in the packings. Let us point 
out here that by long-range order we mean an order that 
extends to the system size when this size tends towards 
infinity [55j . 

In a granular packing, textural order may take various 
forms [5fjj : translational arrangements of grains that 
combine to form patterns, preferential orientation of the 
contact network, preferential orientation of non-spherical 
grains themselves. Each of these aspects is addressed in 
the following paragraphs. 

First, translational arrangements are studied by means 
of the pair correlation function g(r) This function is 
depicted in Fig. [9] for MSP (a) and MPP packings (b) 
and for several values of W. For both packings, local 
order extends over a few particle diameters, slightly less 
for MPP than for MSP due to the higher angularity of 
the former resulting in a loss of rotational symmetry. As 
a consequence, very confined packings exhibit ordering 
over their full size. However, for lower confinement (e.g. 
when W I dh = 10 or more), g(r) tends towards 1 when r 
increases beyond 3 eft, indicating the absence of transla- 
tional long range ordering within our packings. 

Next, the existence of preferential orientations of the 
contact network is investigated. For this purpose, Fig.fTUl 
displays for various confinements 2D representations of 
the distributions of contact orientations in MSP away 
from the bottom plane and the free surface. Given 
the sidewall-induced-layering evidenced in the /i-thick 
boundary layers (see section IIVB[) . contact from the 
boundary layers (Fig. [TUh ) have been dealt with sepa- 
rately from those located in the bulk region (Fig. ITUb ) . 
Note that no bulk region is present in monodisperse pack- 
ings where W/cIl = 5 and, conversely, no boundary layer 
occurs in PBC packings. Furthermore, inside the bound- 
ary layers (see Fig. ITUa). an anisotropy of contact ori- 
entations is visible regardless of the confinement in the 
x, y and z directions, as well as at roughly 60° to the 
x direction in the xy plane and 30° to the y direction 
in the yz plane (corresponding to compact clusters of 
three spheres close to the sidcwalls). This anisotropy is 
fully consistent with the vertical layering of monodisperse 
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FIG. 9. (Color online) Pair correlation functions of MSP (a) 
and MPP (b) for several values of W/cIl- These exhibit local 
order that is stronger for MSP compared to MPP. When r/dh 
is large enough g(r) tends towards one, indicating the absence 
of long range translational order. 



Hence, the anisotropy along the axes x, y and z is miti- 
gated, whilst other preferential orientations correspond- 
ing to various patterns made of small and large grains 
each in contact with the others are generated. Eventu- 
ally, in order to assess the effect of grain shape on the 
existence of preferential orientations of the contact net- 
work, 2D representations of the distributions of contact 
orientations in MPP away from the bottom plane and 
the free surface are represented in Fig. [12] for various 
confinements. As for sphere packings, contacts located 
in boundary layers have been dealt with separately from 
those located in the bulk region. Inside the boundary 
layers, an anisotropy of contact orientations is visible re- 
gardless of the confinement in the x, y and z directions 
(sec Fig. [T2"a ) , and this anisotropy may be explained by 
the wall-induced layering just like for sphere packings. 
In the bulk region, Fig. [T2b shows that pinacoid pack- 
ings exhibit isotropic contact orientation distributions in 
the xy plane, but not along the z axis where, unlike for 
sphere packings, some anisotropy is visible even for pack- 
ings with PBC. This anisotropy may be explained by the 
deposition under gravity protocol, with pinacoids rotat- 
ing around their center of inertia under stcric hindrance 
contraints in order to minimize their potential energy. 

The anisotropy observed in Fig. HTb and WZb respec- 
tively for BSP with x s = 40% and MPP with PBC may 
be calculated and compared to that of the isotropic ref- 
erence state depicted in Fig. [TOb with PBC. Given the 
rotational symmetry of the contact normal orientation 
distributions in the xy plane (see Fig.llOb.lllb and 1 12b). 
the anisotropy may be quantified using a second order de- 
velopment of the contact orientation probability density 
function P(n) (see [43| for details): 



a(3cos 2 6»- 1)], 



(9) 



where: 



packings close to the sidewalls, with a larger peak in the 
y direction due to the high proportion of sidcwall-sphere 
contacts. Unsurprisingly, in the bulk region, Fig. 110b 
shows that the distribution of contact orientations re- 
mains isotropic for all these confinements. In order to 
assess the effect of polydispersity on the existence of pref- 
erential orientations of the contact network, 2D repre- 
sentations of the distributions of contact orientations are 
depicted in Fig. [TT] for W = 5c?L-thick packings, away 
from bottom plane and free surface. Observe that the 
substitution of xs = 40% by mass of small particles for 
large ones results in the emergence of a central 2di-thick 
quasi-isotropic bulk region (see Fig. II lb ) . Furthermore, 
note that the boundary layers remain anisotropic (see 
Fig. \TTh ) although the presence of small particles in be- 
tween large ones tends to disturb the vertical layering of 
the latter (because the centers of inertia of small particles 
are not necessarily coplanar with those of large particles). 



• a = 5(F 3 — i*\)/2 denotes the branch vector coeffi- 
cient of anisotropy derived from eigenvalues F3 and 
Ft of the fabric tensor [U [H HI] 

• 9 denotes the polar coordinate in the xy plane 

This coefficient may vary from (perfectly isotropic pack- 
ing) to 2.5 (perfectly anisotropic packing). Table IVlgath- 
ers values of the branch vector coefficient of anisotropy 
calculated for MSP with PBC, in the bulk region of BSP 
with x s = 40% and for MPP with PBC. These values 
show no significant differences between the bulk region 
of BSP with x s = 40% and MSP with PBC. Further- 
more, the coefficient of anisotropy of branch vectors ob- 
tained for MPP with PBC remains below 0.1, denoting a 
rather small anisotropy. Last, observe in Fig. ll3l showing 
the distributions of orientations of simple (face- vertex), 
double (face-edge) and triple (face-face) contacts that 
the vertical anisotropy in MPP with PBC is identical 
no matter the contact type. As a consequence, no long- 
range contact orientation anisotropy is generated in our 
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FIG. 10. (Color online) Contact normal orientation distribu- 
tions for MSP inside the boundary layers (a), and inside the 
bulk region (b). Upper half of each chart (e.g. from to 180°) 
corresponds to the xy plane, while lower half corresponds to 
the yz plane. Several gap widths are considered. 



frictionlcss grain packings by the grain deposition proto- 
col used, except a weak anisotropy generated in pinacoid 
packings along the z axis. 

Finally, the orientations of particles that are not sym- 
metric by rotation may also be a source of anisotropy 
within the packing. To detect a preferential orientation 
of pinacoids in such a packing, one may use the nematic 



FIG. 11. (Color online) Contact normal orientation distri- 
butions for MSP (x s = 0) and BSP (xs = 40%) inside the 
boundary layers/closer than 1.5di to a sidewall (a), and in- 
side the bulk region (b). Upper half of each chart (e.g. from 
to 180°) corresponds to the xy plane, while lower half cor- 
responds to the yz plane. Gap width is W — 5rf_L and bulk 
region coincides with particles located at least 1.5d_L away 
from sidewalls. 



order parameter Qoo- Here, we recall briefly how this pa- 
rameter can be determined (for details, refer to (57l.[5a|). 
For each particle, if we call ~it ,1? ,~v$ the unit vectors of 
its base of inertia (which, in our case, align with its axes 
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FIG. 12. (Color online) Contact normal orientation distribu- 
tions for MPP inside the boundary layers (a), and inside the 
bulk region (b). Upper half of each chart (e.g. from to 180°) 
corresponds to the xy plane, while lower half corresponds to 
the yz plane. Several gap widths are considered. 



TABLE V. Branch vector coefficient of anisotropy calculated 
for MSP with PBC {W/d L = 20), in the bulk region of BSP 
with x s = 40% (W/d L = 5) and for MPP with PBC {W/d L = 
20). 

Configuration MSP (PBC) BSP (W/d L = 5) MPP(PBC) 
a 0.028 ±0.011 0.032 ±0.005 0.081 ± 0.032 



135^ 
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FIG. 13. (Color online) Distribution of the orientations of 
simple (a), double (b) and triple (c) contacts for MPP in the 
yz plane. Periodic Boundary Conditions are used in the x 
and y direction. 



of symmetry) we can define the following tensor [56 
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with a, = 1,2, 3 



and where n is the number of particles and S the Kro- 
necker symbol. We apply the same definition with Q^L 
and Q™p- From those tensors, the nematic order param- 
eter can be determined: 

Qoo = -Q z " -Z ', 

where Z is the eigenvector corresponding to the larger 
eigenvalue of the three tensors Q uu , Q vv and Q ww . Q zz is 
the corresponding diagonalizcd tensor. By construction, 
this parameter varies between and 1. For each of our 
pinacoid packings, two values of the nematic order pa- 
rameter have been calculated, one corresponding to par- 
ticles located in a boundary layer, and one correspond- 
ing to particles located in the bulk region. Whatever 
the confinement, the nematic order parameter ranges be- 
tween 0.05 and 0.07 in the bulk region, which is quite 
low and shows the absence of privileged grain orienta- 
tion, whereas it is slightly higher in the boundary layers 
(between 0.105 and 0.128). 

As a conclusion of this subsection, frictionless grain 
packings used in the present work do not exhibit sig- 
nificant long range order, except a weak anisotropy of 
the contact orientation distributions observed in pina- 
coid packings along the z axis. Furthermore, sidewalls 
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induce order close to their location that, in the case of 
very confined packing, propagates over the whole system. 



B. Coordination number 

Figure [14] shows the variations of the coordination 
number with d L /W for MSP, BSP and MPP. Each value 
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FIG. 14. (Color on-line) Coordination number as a function 
of <1l/W for MSP and MPP. Error bars on pinacoid pack- 
ings results denote the standard deviation (not represented 
for sphere packings because errors are smaller than symbol 
size). The linear relationship between Z and dhjW suggests 
that the geometrical model, initially derived for the packing 
fraction, is also valid for the coordination number. 



is averaged over three simulations and the error bars 
denote the corresponding standard deviation. Prelimi- 
nary examination of our results obtained with bi-pcriodic 
boundary conditions (unconfined state with d^/W — > 0) 
suggests the following remarks: for sphere packings, the 
calculated coordination number is 6.027 ±0.012, which is 
very close to the 6.073 ± 0.004 value calculated by [48| in 
the RCP state. For pinacoid packings, the calculated co- 
ordination number is 8.581 ±0.068. Although no study of 
pinacoid packings could be found in the literature, such 
a high coordination number value has already been ob- 
served in disordered packings of particles having a similar 
shape (8.6 ±0.1 calculated by |59| for packings of tetrae- 
dra upon extrapolation to the jamming point). 
When confinement increases, the coordination number 
decreases linearly for both MSP and MPP, which is con- 
sistent with the linear decrease of the solid fraction ev- 
idenced in Fig. [To] Though, MPP coordination number 
values tend to be more scattered than MSP ones, which 
could be due to a combination of finite packing size ef- 
fects with the higher level of interpenetration observed in 
pinacoid packings. Nevertheless, the aforementioned lin- 
ear relation between Z and 1/W suggests a generalization 
of the geometrical model to the coordination number. 



For this purpose, let us define Zb u ik and Zbl, respec- 
tively the coordination number for the bulk region and 
the coordination number for the boundary layers. By 
writing the coordination number as the average of Z^uik 
and Zbl weighted by the thicknesses of their respective 
zones (resp. W — 2hz and 2hz) we obtain 



Z = 



W - 2/7. ; 



W 



Jbulk 



2h z 
W 



with Cz = 2hz(Zbuik - Zbl)- 

Figure [14] also shows that the influence of polydis- 
persity on packing coordination number Z decreases to 
zero when the confinement diminishes, which is consis- 
tent with [oOj . Indeed, in the unconfined state, the lack 
of contacts of small spheres with others (due to the steric 
hindrance of large ones) is compensated by the excess of 
contacts of large spheres with small ones. 

To investigate the coordination number decrease with 
increasing confinement, Fig. [T5] depicts coordination 
number profiles in the y direction (normal to the side- 
walls) for sphere and for pinacoid packings. Each of 
these profiles is averaged over three simulations and is 
determined upon subdividing the packing into slices per- 
pendicular to the y direction and calculating for each 
slice the average number of contacts per particle having 
its center of inertia in the slice. In confined state, all 
these profiles evidence a central zone where the coordi- 
nation number is almost unchanged compared to the un- 
confined reference state (except for sphere packings with 
W = 5di), and two "drop zones" in contact with the 
sidewalls where the coordination number symmetrically 
drops by 1.4 (for sphere packings) to 1.5 contacts (for 
pinacoid packings) from their respective unconfined ref- 
erence state. The thicknesses of these drop zones look 
identical to that of the boundary layers described in the 
geometrical model [l8|, leading to the same conclusion 
that grain angularity mitigates the effect of sidewalls on 
the coordination number drop in their vicinity. 
To confirm this observation, we may consider the geo- 
metrical model and compare = C/<f>buik with £z = 
C z /Z bulk . For MSP, we obtain Q = 0.304 and (z = 
0.305 and for MPP g, = 0.186 and ( z = 0.193. Note that 
the values of ^ and C,z are also comparable in the case of 
BSP. The strong correlation between those two quantities 
shows that the propagation of the confinement effect is 
comparable for the two studied quantities: <\> and Z . As 
described in subsection III D[ pinacoid packings incorpo- 
rate simple, double and triple contacts and it is of interest 
to investigate the effect of confinement on their respective 
distribution. Therefore, Fig. [If)] depicts the coordination 
number profiles of MPP for simple (a), double (b) and 
triple contacts (c). Like in Fig. Qj)] all these profiles ev- 
idence a central zone where coordination number values 
are almost unchanged compared to the unconfined refer- 
ence state (except for sphere packings with W = 5di). 
These values are Z s ss 5, Zd ~ 1.7 and Z t ~ 1.8 respec- 
tively for simple, double and triple contacts. 

In order to check the relevance of these coordination 
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FIG. 15. (Color online) Coordination number profiles (along 
y) for MSP (a) and for MPP (b) for several gap width. These 
profiles evidence a constant central zone and two drop zones 
in contact with the side walls. 



number values, one shall observe that packings of n fric- 
tionless rigid grains at equilibrium obey the following re- 
lation [49[ between the degree of hypostaticity fco, the 
degree of hyperstaticity ho, the number of contacts that 
carry forces N c = \-{Z s + 2.Z c i + 3.Z t ), and the number 
of degrees of freedom Nf = 6n of the packing: 



2.5 



12 + 2.^ = Z, 
n 



N f + h = N c 
- 2.Z d + 3.Z t + 



2> (11) 

n 



If we assume that the pinacoids in our packings are ran- 
domly oriented, which seems reasonable according to the 
values of the nematic order parameter (see section IV^f . 
then no motion is possible without generating work in 
the contacts network, which means that the degree of 
indeterminacy of contact forces in the packing is zero, 
therefore ko should be set to in equation [TT] Upon 
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FIG. 16. (Color online) Coordination number profiles (along 
y) of MPP for simple (a), double (b) and triple contacts (c) 
for confined (W B [5<f.L, 20g(l]) and unconfined packings. 
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incorporating in equation [TT] the aforementioned coordi- 
nation number values as well as that of fcg, one obtains: 



ho 



13.. 



12 = 1. 



(12) 



Observe that, for isostatic pinacoid packings (e.g. ho = 
0), equation [11] would lead to Z s + 2.Z d + 3.Z t = 12. 
Here, it is clear that Z s + 2.Z d + 3.Z t > 12. The le vel 
of interpenetretation calculated in subsection MI Bl to- 
gether with the finite size of packings may lead to a sum 
Z s + 2.Z d + 3.Z t slightly higher that 12, but it is doubt- 
ful that this sum would reach 13.8 upon this sole expla- 
nation. The presence of hyperstaticity in our pinacoid 
packings seems more realistic and at least consistent with 
equation [12] and with our finding of as much as 2 triple 
contacts per grain (Zt ~ 1.8). Although interesting, fur- 
ther investigation of the presence of hyperstaticity falls 
beyond the scope of the present paper. 

Coming back to Fig. [TH] the profiles show that MPP 
exhibit more simple contacts than the sum of double and 
triple contacts. They also evidence that confinement pri- 
marily impacts the simple contact profiles, whereas dou- 
ble and triple contact profiles remain unchanged. As 
a consequence, the vicinity of sidcwalls is not a privi- 
leged location for edge-to-face or face-to-face contacts, 
although a drop in the simple contact profiles tends to 
make them look over-represented. 

For a fixed confinement, Fig. [171 shows the influence of 
polydispersity on coordination number profiles in the y 
direction (normal to the sidewalls) for sphere packings. 
As before, each of these profiles is averaged over three 
simulations and is determined upon slicing the packing 
perpendicular to the y direction. As observed in Fig. [TJd, 
an increasing polydispersity does not seem to impact 
the bulk region, but rather reduces the thickness of the 
boundary layers, hence mitigates the effect of sidewalls 
confinement on the coordination number. 

Finally, like for the solid fraction, we have not ob- 
served any effect of the grain- wall friction coefficient on 
the coordination number of MSP. 



VI. CONCLUSION 



short-range order in their vicinity. Except a weak con- 
tact orientation anisotropy observed with pinacoid pack- 
ings in the vertical direction, no long-range order was ob- 
served in our packings. We have demonstrated that both 
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FIG. 17. (Color online) Coordination number profiles (with 
y) for BSP. The gap between sidewalls is W — 5<1l- 



the polydispersity and the angularity of grains lower the 
confinement effect. This effect has been observed for the 
solid fraction and for the coordination number. Our re- 
sults have shown that the geometrical model (l8l - l2l| that 
captures the linear evolution of the solid fraction versus 
—1/W is valid for sphere packings as well as for pinacoid 
packings and that it holds whatever the packing polydis- 
persity 

Interestingly, this model, initially derived for the packing 
fraction can be extended to capture the effect of confine- 
ment on the coordination number. The characteristic 
length quantifying the effect of the sidewalls is found to 
be the same for those two quantities. 
Finally, we have shown that the effect of wall friction is 
negligible, indicating that the major influence of the con- 
fining sidcwalls is geometric. 

Several perspectives arise from this study, among which 
the need to address with more details the presence of 
hyperstaticity in our packings of frictionless pinacoids. 



In this work, we have shown how a confining boundary 
alters the solid fraction as well as the internal structure 
of static frictionless granular materials compacted under 
their own weight using the non-smoothcontact dynam- 
ics simulation method. We did not restrict ourselves to 
sphere packings but extended our work to packings made 
of a particular type of polyhedra: pinacoids. 
As previously reported, the presence of sidewalls induces 
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